A projection operator approach to the Bose-Hubbard model 
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We develop a projection operator formalism for studying both the zero temperature equilibrium 
phase diagram and the non-equilibrium dynamics of the Bose-Hubbard model. Our work, which 
constitutes an extension of Phys. Rev. Lett. 106, 095702 (2011), shows that the method provides 
an accurate description of the equilibrium zero temperature phase diagram of the Bose-Hubbard 
model for several lattices in two- and three-dimensions (2D and 3D). We show that the accuracy of 
this method increases with the coordination number zo of the lattice and reaches to within 0.5% of 
quantum Monte Carlo data for lattices with zq = 6. We compute the excitation spectra of the bosons 
using this method in the Mott and the superfluid phases and compare our results with mean-field 
theory. We also show that the same method may be used to analyze the non-equilibrium dynamics 
of the model both in the Mott phase and near the superfluid-insulator quantum critical point where 
the hopping amplitude J and the on-site interaction U satisfy zqJ/U <C 1. In particular, we study 
the non-equilibrium dynamics of the model both subsequent to a sudden quench of the hopping 
amplitude J and during a ramp from J; to J/ characterized by a ramp time r and exponent a: 
J(t) — Ji + (Jf — Ji)(t/T) a . We compute the wavefunction overlap F, the residual energy Q, the 
superfluid order parameter A(t), the equal-time order parameter correlation function C(t), and the 
defect formation probability P for the above-mentioned protocols and provide a comparison of our 
results to their mean-field counterparts. We find that Q, F, and P do not exhibit the expected 
universal scaling. We explain this absence of universality and show that our results for linear ramps 
compare well with the recent experimental observations. 

PACS numbers: 03.75.Lm, 05.30.Jp, 05.30.Rt 
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I. INTRODUCTION 

Ultracold bosonic atoms in optical lattices provide 
us with an unique setup to study properties of bosons 
near a Mott insulator-superfluid (MI-SF) quantum crit- 
ical pointii^. A careful analysis of such experimen- 
tal bosonic systems in optical lattices show that their 
low-energy properties are well described by the Bose- 
Hubbard model 3 , which has already been theoretically 
studied using both analytical^— and numerical^ tech- 
niques. The presence of such an experimental test bed 
has led to a plethora of new theoretical studies on 
the modeler—. Many of the earlier analytical studies 
have concentrated on obtaining the phase diagram of 
the model by using mean-field theory^, excitation en- 
ergy computation 6 , and strong-coupling expansion for 
the boson Green function^. The results obtained by 
these methods have been compared to extensive quan- 
tum Monte Carlo (QMC) data^ii. Out of these meth- 
ods, the strong-coupling expansion^ (excitation energy 
computation 6 -) and the NPRG approach^ provide the 
closest match to QMC data in 2D (3D). 

Recently, it has been realized that such ultracold 
bosonic systems also allow us easy access to the non- 
equilibrium dynamics of its constituent atoms near the 
MI-SF quantum critical point. The theoretical study 
of such quantum dynamics on various models has seen 
great progress in recent years^. Most of these works 
have either restricted themselves to the physics of inte- 
grable and/or one-dimensional (ID) models or concen- 
trated on generic scaling behavior of physical observable 



for sudden or slow dynamics through a quantum critical 
pointi^—. However quantum dynamics of specific exper- 
imentally realizable non-integrable models in higher spa- 
tial dimensions and strong coupling regime has not been 
studied extensively mainly due to the difficulty in han- 
dling quantum dynamics of plethora of states in the sys- 
tem's Hilbert space. The Bose-Hubbard model with on- 
site interaction strength U and nearest neighbor hopping 
amplitude J, which provides an accurate description for 
ultracold bosons in an optical lattice, constitutes an ex- 
ample of such models. Most of the studies on dynamics of 
this model have concentrated on d = 1—, weak coupling 
regime 2 ^, and mean-field order parameter dynamics fol- 
lowing a sudden ramp in the strong coupling regime^ir— . 
Recent experiments 2 clearly necessitate computation of 
dynamical evolution of several other quantities in higher 
dimensional Bose-Hubbard model in the strong-coupling 
regime (U ^> J) beyond the mean-field theory and for 
arbitrary ramp time r. However, none of the works men- 
tioned above presents an analysis of the non-equilibrium 
dynamics of the model beyond mean-field theory. 

More recently, the authors of Ref. H3 have developed 
a theoretical formalism which enables one to analyze the 
dynamics of the Bose-Hubbard model beyond mean-field 
theory near the MI-SF critical point^l The method 
uses a projection operator technique which enables us 
to account for the quantum fluctuations over the mean- 
field theory perturbatively in Jf/U(J(t)/U) and there- 
fore yields accurate results as long Jf/U(J(t)/U) <C 1 
for sudden(ramp) dynamics. This allows one to treat 
sudden and slow ramps at equal footing near the MI- 



SF quantum critical point. As shown in Ref. [U, the 
projection operator method yields an accurate phase di- 
agram and also provides an estimate of dynamically gen- 
erated defect density which shows a qualitatively rea- 
sonable match with recent experimental results^. In the 
present work, we extend these results in several ways. 
First, we present a generic phase diagram of the Bose- 
Hubbard model as a function of the lattice coordination 
number zq and compare these results to the available 
QMC data for several one-, two-, and three-dimensional 
lattices. Our comparison demonstrates that the accu- 
racy of the projection operator technique increases with 
zq reaching about 0.5% of the QMC data for lattices with 
zo = 6. Second, we compute the excitation spectrum us- 
ing our approach and show that it yields the gapless Bo- 
goliubov and gapped amplitude modes in the SF phase 
and the gapped particle and hole excitation modes in the 
MI phase. Third, we study the dynamics of the model 
for non-linear ramp of the hopping parameter J from Jj 
to J f characterized by a ramp time r and exponent a: 
J(t) — Ji + (Jf — Ji)(t/r) a . We compute the fidelity sus- 
ceptibility F, nearest-neighbor correlation between the 
bosons B, the defect formation probability P, and the 
residual energy Q of the system following such a proto- 
col and show that our result reproduce those of Ref. [H 
for a = 1 as a special case. We also find the value of 
the optimal a which leads to minimal defect production 
for fast quenches (small t). Finally, we also compute 
the order parameter A(t), the order-parameter correla- 
tion function C(t), the wavefunction overlap F, and the 
residual energy Q subsequent to a sudden quench, discuss 
their properties, and provide explicit analytical expres- 
sions for A(t) and Q. We also provide a detailed com- 
parison of the behavior of A(t) with that obtained from 
Gutzwiller mean-field theory. 

The plan of the rest of the work is the following. In 
Sec. [Ill we develop the projection operator formalism and 
apply it to obtain the equilibrium phase diagram of the 
Bose-Hubbard model for arbitrary zq and compute its 
excitation spectrum. This is followed by Sec. IIII1 where 
we discuss the dynamics of the model both for sudden 
quench and non-adiabatic ramp of the hopping amplitude 
J. Finally we discuss our results and conclude in Sec. lIVI 



II. FORMALISM AND EQUILIBRIUM PHASE 
DIAGRAM 



In this section, we provide a detailed exposition of the 
projection operator formalism. In Sec. Ill A[ we compute 
the MI-SF phase boundary using this formalism for vari- 
ous lattices while in Sec. Ill Bl we compute the low-energy 
excitation spectra of the MI and the SF phases. 
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FIG. 1: (Color online) Schematic representation of the Mott 
state with n — 1. (b) Typical hopping process mediated via 
T/. (c) Hopping process mediated via T®. Notice that the 
states in (c) become part of the low-energy manifold near the 
critical point, while those in the right side of (b) do not and 
are always at an energy U above the Mott state. 

A. Phase boundary 

The Hamiltonian of the Bose-Hubbard model is 

H = T+H , T=J2 -Jb%> 
<«"} 

H = ^2[-fin r + ^n r (hr - 1)] (1) 

r 

where b r (n r ) is the boson annihilation (number) operator 
living on site r of a d-dimensional lattice with coordina- 
tion number z$ = S( r ') an< ^ the chemical potential /i 
fixes the total number of particles. The exact solution of 
H is difficult even numerically due to the infinite dimen- 
sionality of the Hilbert space. A typical practice is to 
use the Gutzwiller ansatz \ip) — J7 r \ n ) > an< ^ s °l ve 

(r) 

for Cn keeping a finite number of states n around the 
Mott occupation number n — n. This yields the stan- 
dard mean-field results with c n = c n for homogeneous 
phases of the models. 

To build in fluctuations over such mean-field theory, 
we use a projection operator technique^. The key idea 
behind this approach is to introduce a projection opera- 
tor 

Pi = \n)(n\ r x \n)(n\ r > (2) 

which lives on the link t between the two neighboring 
sites r and r' of the lattice. The hopping term T can 
then be formally written as 

T = E^ = E( T °+^) 

T$ = PtT t Pt, Tj = (P i T i + T i Pt), (3) 
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where P^~ = (1 — Pi). The advantage of the decomposi- 
tion given by Eq. [3] is that it distinguishes between low- 
and high-energy tunneling processes as shown schemati- 
cally in Fig. Q] for n — 1 . The existence of a low-energy 
subspace for the model becomes more evident by rewrit- 
ing the Bose-Hubbard Hamiltonian in the more conve- 
nient way H — Ho + Hi, where 



n 



(4) 



We can then define the low-energy subspace to be a set 
of states which are separated from the ground state of 
Ho by energies O(J). These set of states can not be 
connected to each other by Hi- For any two members, 
\ni) and |n 2 ) of this set, one has {n\\Hi\n2) = 0. In 
other words, Hi acting on any state \ni) in these low- 
energy subspace yields a state \n'i) which is necessarily 
separated from ground state of Ho by an energy 0(U). 
Note that the states which are member of the low-energy 
subspace depend on the value of J/U. For example, the 
states schematically represented in panel (c) of Fig. Q] 
become members of the low-energy subspace near the 
MFSF quantum critical point where J ~ J c ; however, 
these states do not belong to the low-energy subspace for 
J = 0. 

In what follows, we shall use the projection operator 
technique to systematically chart out the effective low- 
energy Hamiltonian by eliminating Hi from H to O(J). 
The canonical transformation operator S which achieves 
this can be written as 



S= S[J] =J2i[Pt,T e ]/U. 



(5) 



It can be easily checked that [iS, Ho] = — J2e so * na * 
the transformation eliminates Tg[J] up to first order in 
zqJ/U. A standard expansion in zqJ/U then leads to 
the effective Hamiltonian H* = exp(iS)H exp(-iS) to 



0(zp 2 /U) 



H* 



Ho 



PtT f Pi 



PeTfPe - T e P e T e 



\P t Tf+TfP t 



\T f P t P„T,j - P,T,T,,P,j 



PtTfTt, 



h.c. 



T t P t T t 



(6) 



Note that the second order terms in H* involves effec- 
tive hopping processes between adjacent links leading to 
spatial correlation between next-nearest neighbor sites; 
higher order terms in zqJ/U systematically build such 
correlations between further neighbors. In this work, we 
restrict ourselves to 0[(zoJ/U) 2 ]. 

Using H* one can now compute the variational ground 
state energy 

E = ^\H\f)-W\H*W) + 0{zlJ a /U% (7) 
where \ip'} = exp(iS)\il>}, and we use a Gutzwiller ansatz 

\ip') = Yl T J2 n fn \n), for the variational wavefunction 
|V>) in the Mott limit (S, J = 0), where \tp) = \ip'). Note 
that \ip) is not of the Gutzwiller form; it incorporates 
spatial correlation via exp(zS') factor. To obtain the vari- 
ational energy E in terms of the coefficients, we define 
the fields 



n n 
n 

£y (n + l) (n + 2)/*M/« 2 (8) 



Using the expressions of (p T and $ r in Eq.[8l one obtains, 
after some algebra, the expression for the variational en- 
ergy E = E[{f n }; J] to be 



E = ]T[-/in + Un(n 1 )/2] | | 2 - J £ {^ r , - 2^-iVr'n + M(n + l)/tf [|/«| 2 |jf >| 2 



<rr'> 



-l/Srn/ni 2 -/^/^!/^?/^] +2J/Um:. n ^ rlfl }-J 2 /U £ {2^[< rl _ 1 (n+l)|/r ) | 2 

(rr'r") 



-<Pr,n-l(n + ±)\fn + \\ 2 <Pr" + Vrn n 
+2»^rt*r' 1 fi-lVj».fi-ll, 



\A-\\ 2 - \ffPl 2 } Wr«« + Wr.n-xCft + 1) [|/£i| 2 - |/f >| 2 



(9) 



Note that the first three terms in the first line of Eq. [S] are corrections due to quantum fluctuations. Thus the 
represent the mean-field energy functional, while the rest projection operator method involves a systematic way of 
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incorporating quantum fluctuations over the mean-field 
theory and we expect the results from this method to be 
accurate for larger z$ where mean-field theory provides 
an accurate starting point. 

The Ml-SF phase diagram can be obtained by min- 
imizing E[{f n }\ J] with respect to {/„} or by solving 
ihdt\ip') — H*[J]\ip') in imaginary time2&. In this work, 
we are going to use the former technique and restrict 
ourselves to n = 1. Such phase diagrams for 2D and 3D 
square lattice are shown in Fig. [2fa) and Fig. EJb) re- 
spectively. We note that the match with QMC data^i is 
nearly perfect for 3D (Figj2fb)) where mean-field theory 
provides an accurate starting point. While in 3D the ac- 
curacy with QMC at the tip of the Mott lobe is ~ 0.5%, 
in 2D we find J c /U = 0.055 compared to the QMC value 
0.061-i (red line in Fig.^a)). Here the match with QMC 
is not as accurate as in 3D; however it compares favorably 
to other analytical methods^. 

To provide an accurate comparison of our method with 
other lattices, we note that the nature of the lattice af- 
fects J c only through zq. We also note from Fig.^b) that 
the deviation of J c from QMC value is maximal at the 
tip of the Mott lobe. Thus to elucidate the zq and hence 
the lattice dependence of J c computed by the present 
method, we plot J* lp (i.e., the value of J c at the tip of 
the Mott lobe) as a function of zq in Fig. [3J The compari- 
son of corresponding QMC data for various lattices show 
that the method indeed becomes more accurate with in- 
creasing Zq. 

Before ending this section, we would like to note 
that the inclusion of fluctuation in our method be- 
comes apparent on computing the expectation (T e ) — 
— J{xj}\b\b r i -f-h.cIV') in the MI phase. The mean-field the- 
ory provide a zero result for (Tf) , while the projection 
operator method yields 

(VIW) = (ip'\eMiS[J])TteM-iS[JW) 

= W) - i(V/| [PtTj + T 2 P e - 2T e P e T e ] 

+ [PiTiT e < - T t P t T' t + h.c] (10) 

where we have kept terms up to 0(J 2 /U 2 ) and (£') de- 
notes nearest neighbor links to I. In the MI phase, the 
first term of Eq. [TU1 which is also the mean-field result, 
vanishes, while the second fluctuation contribution from 
the remaining terms yields (T/) = 2J 2 n(n + 1)/U in the 
homogeneous limit. We note that this agrees with fluc- 
tuation calculations of Ref. H 



B. Excitation Spectrum 

To obtain the low-energy excitations, we consider a 
variational form for which corresponds to perturba- 




FIG. 2: (Color online) Phase diagram of the Bose-Hubbard 
model in 2D (a) and 3D (b). The blue dots and blue solid lines 
(black dashed line) indicate the phase diagram obtained by 
the projection operator (mean-field) method. The red squares 
indicate QMC data. 
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FIG. 3: (Color online) The plot of J* lp as a function of zq 
as shown by blue circles. The red square, brown hexagon, 
green triangle, and the blue inverted triangle represents QMC 
data for ID Bose Hubbard model (zo = 2), 2D square lattice 
(zo = 4), 3D cubic lattice (z.o = 6), and 2D triangular lattice 
(zo = 6)respectively. 

tion over the ground state value. This is given by 

w) = nE/n r) w 

r 7i 

fi r) (t) = If^ + Sf^iW^ (11) 

where 6fn\t) represents small perturbation over the 
ground state value /„ and can be expressed in momen- 
tum space as 

<y/W(t) = ule 1 ^-^ +v£e- i0 <- r - ut \ (12) 

Substituting Eqs.[TT1and[T2lin the Schrodinger equation 
ihdt\ip') = H*\ip'), we obtain a set of equations for it£ 
and v£ which is given by 

-*)(£)■ (13) 

Here the Uk and Vk are vectors with components 
and v£ (n = 0, 1, . . . ), respectively and and are 



5 



square matrices with elements A£ n and B™ 1 . Since 
< n < oo, in principle, Eg. 1131 represents an infinite- 
dimensional matrix equation; however, in the strong 
coupling regime where states with n > 3 bosons are 
energetically costly, it is possible to truncate the ar- 
rays Mk and Uk by putting u^,v^ = for n > 3. In 
this case the column vector (itk)^) can be written as 
(lf k , it 2 ., "ki w k ' w k' u k' w k) T - Thus the solution for the 
excitation spectrum reduces to the solution of a 8 x 8 ma- 
trix for each k. In what follows we provide an analytical 
solution for in the MI phase and a numerical plot of 
the excitation spectrum in the SF phase, where the al- 
gebra, for reasons mentioned below, turns out to be too 
complicated to allow a straightforward analytical result. 
We note here that our analysis amounts to generaliza- 
tion of the work in Ref. [27| which provides the excitation 
spectrum of the Bose-Hubbard model using mean-field 
theory 

In the MI phase, fn^ — Sni and hujo can be shown 
to correspond to the ground state energy of the homoge- 
neous MI state as obtained by putting /} ' = 1 in Eq. |H1 
hujo = — /i — Az J 2 /U. Further, one finds that in the MI 
phase the elements A™ 1 ~ S nm . Using the expression of 



wo, these diagonal elements can be calculated to be 

2J 2 z a 



A ™ = M - Jzo (l - 2 ^ 2 ) - 
4J 2 (1 - 2a; 2 ) z„ 



U 



A 11 

I 22 



u 



Al 2 = U - 2/t - Jzq (l - 2a; 2 ) + A 



(l-2x 2 ) 2 z 



00 



A? = ^+3t/-2/, 



(14) 



Note that in the limit J = 0, these elements corre- 
spond to the on-site excitation energies of the different 
\n) states. The off-diagonal elements are given by 



B™ = - 



3V2J 2 Z() ((1 - 2x 2 f z - 1 



Bi b = - 



AJ 2 (1 - 2.x 2 
U 



U 
) z a 



= B 



21 



(15) 



where x = Y^a=i d sm2 (^a/2)/^ an( i we nave se t the lat- 
tice spacing to unity. Diagonalization of the matrix in 
Eq. [13] leads to the excitation spectra given by 



-Ekl 



3(17 - /*), Eh 



k2 



^{AJ 2 z Q + Up) [Up - AJ 2 (Ax 2 - 3)] 
U 



-8J 2 z 



" a Z ° ((1 - 2x 2 ) 2 z -lj + ( J 2 z ((l - 2a; 2 ) 2 z a + 8) + 12 J (l - 2x 2 ) z ol i + A^ 2 

1/2 



1 - 2a; 2 ) 2 z - 1 ) (3J (2a; 2 - l) z - 2/i) + 2U (3 J (2a; 2 - l) 2 Q - 2/i) + 



+ Jz (2x 2 - l) + U - 2/i 



(16) 



Note that E\a corresponds to the hole branch while _Bk3 
corresponds to the particle branch. The energies of these 
branches differ from their mean-field counterparts in Ref. 
l27l via presence of additional fluctuation contribution 
which manifest then through 0(J 2 /U) terms. The plots 
of these excitation energies as a function of x is shown 
in the left panel of Fig. |4] and matches qualitatively with 
its counterpart in Ref. 1271 

A similar analysis for the superfluid phase can easily 
be carried out using the same algorithm described above. 
In this case, it turns out that the analytical expressions of 
A\t ( which now has off-diagonal terms) and B^ are pro- 
hibitively lengthy. We therefore resort to numerical solu- 
tion of Eq.[T3]for several values of x. The result is shown 
in the right panel of Fig. SJ The qualitative features of 
the plot s are again similar to the mean-field results of 
Ref. l27t however quantitative values of physical quan- 
tities such as the velocity of the Bogoliubov mode, v g , 
differ. The difference with the mean- field result comes, 



again, from the presence of 0(J 2 /U) terms in the effec- 
tive action and hence is small near the critical point. 

III. DYNAMICS 

We now demonstrate that the method elaborated in 
Sec. |TT] with minor modification allows one to address 
the dynamics of the Bose-Hubbard model. To this end, 
we are going to assume a protocol where the hopping 
amplitude J = Jit) changes in time from its initial value 
Ji to some final value Jf. The necessary condition for our 
method to yield accurate result, as we shall demonstrate, 
is ZQ,J(t)/U < 1 at all times. 

We begin with the Schrodinger equation for the time 
dependent Hamiltonian %[J(t)\ which is given by 

ihd t \ip)=n[j(t)]\ip) (17) 

The solution of this equations is difficult due to the in- 
finite dimensionality of the bosonic Hilbert space. How- 
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FIG. 4: (Color online) Left panel: Plot of the excitation spec- 
tra -Eki (dashed blue line), E\a (solid purple line) and _Ek3 
(dotted yellow line) as a function of x in the MI phase for 
fi — 0AU and J = 0.15J C . Right panel: Analogous plots for 
the excitation branches in the superfiuid phase showing the 
gapless Bogoliubov mode and the gapped amplitude modes 
for fi = 0.4(7 and J = 1.3J C . 



ever, one notes that the contribution to the dynam- 
ics of the bosons, as long as zoJ(t)/U <C 1, comes 
from a limited set of states which are members of 
the instantaneous low-energy subspace at any given in- 
stant t. To capture the contribution of the states in 
this low-energy subspace, we make a time-dependent 
transformation \ip') = exp(iS[J(t)])\ijj}, which elimi- 
nates T°[J(i)] up to first order from H[J(t)} at each in- 
stant, and leads to the effective Hamiltonian H*[J(t)] = 
exp(iS[J(t)])H[J(t)]exp(-iS[J(t)]). This yields the 
equation 



(ihd t + dS/dt)\^) = H*[J(t)]\if>') 



(18) 



We note that the additional term dS/dt takes into ac- 
count the possibility of creation of excitations during the 
time evolution with a finite ramp rate r" 1 . The above 
equation yields an accurate description of the ramp with 
H*[J{t)] given by Eq.©for J(t)/U < 1. Note that this 
does not impose a constraint on magnitude of r; it only 
restricts Jf/U and Ji/U to be small. Thus the method 
can treat both "slow" and "fast" ramps at equal foot- 
ing. Substituting = FJ r fn we obtain a set of 
coupled equations for the coefficients {/ n } 



ih 



dfk 



to 



ih dJ r i— „( r ) ~ — 



J2 

u 



Vn(n-l)fi%rj r + y/(n + \){n + 2)f^r 



(19) 



where the fields \n\ a r , j3 r * , rj r , £ r *, da r , and d/3 r * have 
to be calculated self-consistently and are explicitly given 
in the Appendix [5] Notice that the last line of Eq. IT51 

(r) 

couples the time derivative of the coefficient /A with 
the coefficients fj.+ 2 , the coupling being proportional to 
J 2 /U. In particular this is different than the standard 
mean-field equations where the time derivative of the co- 
efficient /ri is at most coupled to the coefficients /A±i 



through J. It is worth noting that Eq. [19] conserve the 
total number of particles for any J(t). In what follows, 
we shall obtain a numerical solution of Eq. [19] to address 
the dynamics of a translationally invariant Bose-Hubbard 
model both for sudden quench and non-linear ramp of 
J(t). 



A. Sudden Quench 

In this section, we are going to address the dynam- 
ics of the bosons after a sudden quench of the hopping 
amplitude from Ji (Mott phase) to J/ (superfiuid phase) 
through the tip of the Mott lobe where the dynamical 
critical exponent z — 1. Our main objective here is 
to compute the time evolution of the order parameter 
A r (i) = (■ip(t)\b r \ip(t)) and the order-parameter correla- 
tion function C r (t) = (^(t)\b r b r \tp(t)) - A?(i). We shall 
also consider sudden quenches which start at the critical 
point (Ji — J c ) and end in the superfiuid phase Jf > J c , 
and compute the resultant residual energy Q and the 
wave function overlap F. 

To this end, we begin by noting that for a sudden 
quench, dJ/dt ~ S(t) and thus the first term on the right 
side of Eq.[l9]does not contribute to the subsequent time- 
evolution of the system for t > 0. The time evolution of 
the order parameter A(i) can then be written in terms 
of {f n (t}} by noting that A(t) = (ip' (t)\b' r \ip' (t)) , where 
b' r = exp(iS[Jf])b r exp(— iS[Jf]). One can then express 
to 



A(t) in terms of fn as 



Ar(t) 



.(t) + J/Uj2n[\f., 

to)r 



to 1 2 



i^tl 2 



+ (- + l)[l/i r) | 2 -/S 1 | 2 
- $ T ,n-l <P* r >n + $rn - $ 



r.n— 2 



i (20) 



Note that the first term in Eq. [20] represents the mean- 
field result while the presence of the other terms indicate 
contribution from the quantum fluctuations from mean- 
field theory. The role of such quantum fluctuations in 
the evolution of A r (t) becomes evident in computing the 
equal-time order parameter correlation function C r (t). 
To compute A r and C r , we consider a spatially homoge- 
neous system and solve the Schrodinger equation fEq. U9[) 

for fn = fn (as guaranteed by translational invariance) 
keeping all states for < n < 5 with n = 1. The resultant 
plot of A r (t) ee A(t) is shown in Fig. [^a) [(d)] for J t = 
and J f /J c = 1.02(J//J C = 3.51). We find that near 
the critical point, A(t) displays oscillations with a single 
characteristic frequency^ while away from the critical 
point (Jf/J c = 3.51), multiple frequencies are involved 
in its dynamics. The time period T (Fig. [jjc)) of these 
oscillations near J c is found, as a consequence of critical 
slowing down, to have a divergence T ~ (jj)-o.35±o.05 
leading to zv = 0.35 ± 0.05 for d = 3^^. Finally, we 
plot C r (i) ee C(t) as a function of t for Jf — 1.02J C in 
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Fig.EJV). We find that \C(t)/A 2 (t)\ may be as large as 
0.5 at the tip of the peaks of A(t), which shows strong 
quantum fluctuations near the QCP. 

To compare our results with the order parameter dy- 
namics obtained from the mean-field theory, we solve the 
equations of motion for the time dependent Gutzwillcr 
coefficients f n (t) within a single-site mean-field theory. 
As shown in Ref. j2lj}23|, the mean-field equation reads 



(a) 



ill C JJ± 



(hi 



200 400 600 800 1000 

tu 



200 400 600 800 1000 
tU 



i).fn = -z J(i) A m{ (t)^/nf n 



-AJnfWVrT+T/, 



n+l 



(21) 



where e„ = — (in+Un(n— 1)/2 is the on-site energy of the 
bosons, A mf (t) = Y, n fn-JnVn, and J(t) = (J^(-t) + 
Jf9(t)) for the sudden quench protocol. For the mean- 
field theory the critical point lies at J™ f = 0.028C/ for 
d = 3. To obtain the order parameter dynamics, we 
obtain the values of f n {t) numerically keeping up to n = 5 
states and compute the order parameter A m f (t) choosing 
the same sudden quench protocol used in the projection 
operator approach (see Fig. [5]). The behavior of A m f (i) 
as a function of time is shown in the left panel Fig. |6] 
A comparison of our results with that of the mean-field 
theory can now be made by comparing Figs. [5] and [6] We 
find that although the qualitative nature of A m f(i) and 
A(t) are similar, the periodicity of the oscillations are 
quite different. Further we note that 



CnAt) = (b 2 r ) - A 2 mi (t) 

= V n ( n ~ !)/n-2/« 



Ai f (t) (22) 



is also expected to show qualitatively similar behavior 
to C(t). Thus we conclude that the subsequent dynam- 
ics of the order parameter following a quantum quench 
near a critical point is qualitatively similar in nature to 
what is found from mean-field theory; however, the pre- 
cise quantitative value of, for example, its period of os- 
cillation, receives significant contribution from quantum 
fluctuations. 

Next, we compute the wavefunction overlap F = 
|(^# c }| 2 = |(V^|e lS [ J /]e- l5 [ J <=]|^>| 2 for sudden quench 
starting at the QCP. Here ^>f(ip c ) denotes the ground 
state wavefunction for J = J/(J C ). The residual energy 
Q = (ip c \H[J f ]\ip c ) - E G [J f ], where E G [J f ] denotes the 
ground state energy at J = Jf as obtained by minimizing 
E in Eq. [51 can also be computed in a similar manner. 



(r) 



Using the fact that for \ip' e ) = e lS ^ J ^ \ip c ), ip r 
we find, in terms of the coefficients fn 

Q = E G [J C ]- E G [J f ]~2J5Jn(n+l 



0, 



E 

<«■'> 



[i/i 1 



r)|2 i/f ) i 2 



I f( T ) |2|f{r')|2 f*WfW j-*( r ') f( r ') 1 /rr 
"IJft+ll \Jn-l\ J ri+lJ fi-lJn-1 J n+l\ I u ■ 



A plot of 1 — F and Q for the homogeneous case, as a 
function of SJ for SJ/J C < 0.2 is shown in Fig. \7\ A 



5 

H 



200 




FIG. 5: (Color online) Plot of A(t) (a), and C(t) (b) as a 
function of tU, for J/ = 1.02J C . (c) The time period T of the 
oscillations of A(t). (d) Same as in (a) for Jf — 3.51J C . We 
have set h = 1 for all plots. 
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FIG. 6: (Color online) Plot of A(t) as a function of Ut (h = 1) 
as computed using mean- field theory for J; = 1.02J c mt (left 
panel) and Jf = 3.51J™ (right panel). 
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FIG. 7: (Color online) Plot of F and Q as a function of the 8 J 
for SJ I J c <C 1. The lines correspond to fits yielding a power 
1 - F(Q) ~ (<5J) ri(r2) with ri ~ 0.89 and r 2 ~ 1.9. 
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numerical fit of these curves yields 1 — F ~ 8 J - 89 and 
Q ~ (5J 190 which disagrees with the universal scaling 
exponents (1 - F ~ 5J diy and Q - <5j( d +*>) expected 
from sudden dynamics across a QCP with z = 1—. In 
the next section, we shall study ramp dynamics across the 
quantum critical point and try to understand the reason 
behind such lack of universality in the dynamics of the 
bosons. 

Before ending this section, we note that it is possi- 



ble to compute the evolution of the correlator Bi — 
{{b\b rl +h.c.)) (I is the link between sites r and r') after 
a sudden quench from J = Ji to Jf where Ji corresponds 
to the MI state and Jf correspond to either the SF phase 
or the MI phase. The mean- field results for such a cor- 
relation would be zero if Jf corresponds to the MI state 
and |A 2 (£)| if it corresponds to a homogeneous SF phase. 
In contrast, using Eq. [TU]with J = Jf, we find that the 
projection operator approach yields, 



Bp 



^{Kn-2^,n + $r,R^, fl _ 2 ) +2n(n+ 1) IfPW^ | 2 

(fi 4-^ ff^ fW* f ( r ') f ( r ')* , f W f (r)* f 0') f ( r ')* ,,(r) |2| f 0')|2 



u 



23? 



<«-i(n + l)|/ri 2 



-y* n $ rl ,n-i - (fr,n-i(n + l)\fn+i\ 2 ]<PP>}, 



(23) 
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FIG. 8: (Color online) Plot of Bt(t) after a sudden quench 
from Ji to Jf as a function of time for Jf = 0.98J C (left 
panel) and Jf = 1.02J C (right panel). See text for details. 



When Jf correspond to the MI phase, since /«(£) 3> 
fn^n(t), we find that Bi shows very small oscillations 
around the base value AJ f\ff L \ 4L /U . In contrast, it dis- 
plays significant oscillation in the SF phase. These be- 
haviors, in the homogeneous limit, are sketched in the 
left and right panels of Fig. [5] for Jf = 0.98J C (left panel) 
and I.02J C (right panel). We note that the behavior of 
Bi in the MI phase is qualitatively different from the 
mean-field result which correspond to the first term in 
the right side of Eq. 



B. Non-linear Ramp 

In this section, we address the dynamics of the bosons 
during a ramp of the hopping amplitude J character- 
ized by a rate r _1 and an exponent a: J(t) = Ji + (Jf — 
Ji)(t/ T ) a ■ Note that the system evolves from at t, = 
to Jf at tf — t; consequently as long as we restrict our- 




rU 



FIG. 9: (Color online) Plot of F = 1 - P as a function rU 
(in units of H = 1) for Ji/U = 0.05 (SF phase) and J f /U = 
0.005 (Mott phase) for a — 1 (blue circles), 2 (red squares), 
3 (yellow diamonds), 4 (green triangles), and 5 (blue inverted 
triangles) showing the plateau-like behavior at large r. 



selves to Ji/U, Jf/U -C 1, we expect the perturbative 
projection method to address the dynamics accurately 
irrespective of the values of r and a. Thus the projec- 
tion operator method enables one to address "slow" and 
"fast" and linear/non-linear ramps at equal footing. We 
note at the outset that our results in this section repro- 
duce those in Ref. [13 as a special case for a = 1. 

To address the dynamics, we use Eq. [19] and solve for 

(r) 

fn = f n for translationally invariant systems. This 
enables us to compute the defect formation probability 

P = 1 - |(^ G ^(t/))| 2 = 1 - F, where \i> G ) (|^(*/))) 
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FIG. 10: (Color online) Plot of Q as a function of r for a — 
1..5. All parameters and symbols are same as those in Fig. [9] 
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FIG. 11: (Color online) Plot of dF/dr for small r (tU < 5) 
as a function of a. 



denotes the final ground state (state after the ramp), 
for a ramp from Ji/U = 0.05 (superfluid phase) to 
J f /U = 0.005 (Mott phase) as a function of r. The 
behavior of F and Q are shown in Figs. [UJ and QJJ for 
various representative values of a. We find that both Q 
and F (and hence P) exhibits a plateau like behavior at 
large r. The slope of both F and Q for small r depend 
on the ramp protocol through the exponent a; however, 
the asymptotic values of these quantities at large r is in- 
dependent of a. The plot of dF/dr for tU < 5 (where F 
is approximately linear in r as can be seen from Fig. [UJ) 
is shown in Fig. [TTJ The slope decreases monotonically 
with a for large a which indicates that F (and similarly 
Q) saturates at larger values of r with increasing a. The 
slope is maximal, indicating minimal initial defect pro- 
duction, for a = 1.5. 

It is clear from the plots that both P and Q do not 
display universal scaling as expected from generic theo- 
ries of slow dynamics of quantum systems near critical 
pointed. This seems to be in qualitative agreement with 
the recent experiments presented in Ref. |2|, where linear 
ramp dynamics of ultracold bosons from superfluid to the 
Mott region has been experimentally studied. Indeed, it 
was found, via direct measurement of parity of n per site, 



that F displays a plateau like behavior similar to Fig. |9j 
Such a lack of universality in the dynamics can be qual- 
itatively understood from absence of contribution of the 
critical (k = 0) modes. In the strong-coupling regime 
(zqJ/U -C 1), the system can access the k = modes 
after time T which can be roughly estimated as the time 
taken by a boson to cover the linear system dimension 
L. For typical small J (U = 1) in the Mott phase and 
near the QCP, T ~ 0(Lh/J) can be very large. Thus 
for t < T ', the dynamics, governed by local physics which 
is well captured by our method, do not display critical 
scaling behavior. We note that our theory which is based 
on building on spatial correlation order by order in pow- 
ers of zoJ/U shall not easily capture the physics associ- 
ated with long-range spatial correlation near the critical 
point and will deviate from experimental results for much 
slower ramp rates. It seems, however, that achieving such 
low ramp rates for the present system in the Mott phase 
can be experimentally challenging. 



IV. DISCUSSION 

In conclusion, we have presented a projection oper- 
ator formalism that describes in a semi-analytical way 
both the phase diagram and non-equilibrium dynamics 
of the Bose-Hubbard model. It produces a phase dia- 
gram which is nearly identical to the QMC results in 
3D, allows for a computation of the low-energy excita- 
tion spectra of the system, and yields semi-analytical in- 
sight for several quantities such as F, Q, A(t), P, and 
C(t) for non-equilibrium dynamics. Its prediction for P 
for a slow ramp matches qualitatively with recent exper- 
iments. The method, in principle, can be generalized to 
any strongly correlated systems which allows perturba- 
tive treatment of fluctuations. We leave such considera- 
tions for future study. We also note that studying finite 
temperature physics of the Bosons with our method also 
poses an interesting theoretical challenge. For now, we 
can only estimate the range of physical temperatures T 
for which the T = theory is accurate. For a typical 
lattice depth in the Mott or critical regime, one can esti- 
mate U ~ 2 kHz ~ 200nKi. This yields, in 3D, a melting 
temperature T* ~ 0.2J7 = 40nK for the Mott phase and 
critical temperature T c ~ zq J c ~ 35nK for the SF phase 
at Mott tip2£. This necessitates T -C T c , T* to be a few 
nano-Kelvins which is well within the current experimen- 
tal limit ~ lnK22. 
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Appendix A: Explicit form of the fields 

In this section we provide explicit expressions for the 
fields used in Eq[TjJ] In what follows we define <5y to be 
the Kronecker delta and Ay = 1 — <5y . The fields that 
multiply the time derivative of the tunneling coefficient 
are given by 

dai = (Snnfbn — $n,n+i<Pb,n-i) i (Al) 
Wi = £ (SnnVU-i ~ K*-l<PU) ■ (A2) 

The phase factor Xn is given by 
U , N 2J 2 



\ : -fJ-n + —n(n - 1) - —S nn n(n + 1) £ |/| o) 



/(<) : 

IT 



(a), <c>, 



.* </>cn 



+ (n + l)(S nn ~ £n,n+l)^>a,n-l 



^n(n + 1) [<5„, ft+1 ]T |if \| 2 + ( 5„, fi _ 1 £ l/W 

<a>, (a)i 



(a) |2' 



Notice that gives only a phase factor in real time 
dynamics and is therefore negligible; however, it is im- 
portant in the imaginary time evolution. The fields a r 
and P*, which couples d t fn^ to f^2i linearly in J are 
given by: 



(A3) 



on 



and 



<<*>i 

— 5n,n+l (fa 



u a 
- 1 - — c« 



(A4) 



J 



g*(t) 



(«>» 



^ a 



^n(< ft -x-^C (i) )]- (A5) 



In contrast, the fields ry r and £ r *, which couple d t fn^ to 
/n±2 to 0(J 2 /C7) are given by 



^ ^ "j^nn^an d" ( ^7i,n+2 < ^ ) a,n — 2 ^n.n+l^djn — 1 



(a) 



^n,n+l A ac ^ 

{c)i 



fan + <Pa,n-l l^c - fanfc,n-l 



}■ 



and 



+ E A ac(<yn,n-2¥'an + <W¥>a,n-l) V* 

(c>< 

(A6) 

where we have introduced the quantities 

A a ] = E A«[(*o,ft-a-*a > fl-l)¥'6ft+(*a«-*a > fl-l)p6,a_i 
(*>« 

+(n+l)(|/^ ) | 2 -|/ii| 2 )^-t+n(|/i o) | 2 -|/t ) 1 | 2 )^], 

(A7) 

B« = £ A b , [(ft + 1) |2 _ 1/(^2) _ 

+ ($an - $a,n-l) ft + $a,n-lftn , (A8) 

and 

CW = ^A«[n(|f | 2 - l/^l 2 )^ - ^) 

W a 

+ ($a,n-2 - <&a,n-X)fl + ^a,n-lft,n-l > ( A9 ) 

for notational convenience. 
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